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Generalized Langevin dynamics (GLD) arise in the modeling of a number of systems, ranging 
from structured fluids that exhibit a viscoelastic mechanical response, to biological systems, 
and other media that exhibit anomalous diffusive phenomena. Molecular dynamics (MD) 
simulations that include GLD in conjunction with external and/or pairwise forces require the 
development of numerical integrators that are efficient, stable, and have known convergence 
{Sj properties. In this article, we derive a family of extended variable integrators for the Gener- 

alized Langevin equation (GLE) with a positive Prony series memory kernel. Using stability 
and error analysis, we identify a superlative choice of parameters and implement the cor- 
responding numerical algorithm in the LAMMPS MD software package. Salient features of 
the algorithm include exact conservation of the first and second moments of the equilibrium 
& velocity distribution in some important cases, stable behavior in the limit of conventional 

Langevin dynamics, and the use of a convolution-free formalism that obviates the need for 
explicit storage of the time history of particle velocities. Capability is demonstrated with 



respect to accuracy in numerous canonical examples, stability in certain limits, and an ex- 
emplary application in which the effect of a harmonic confining potential is mapped onto a 
■^- memory kernel. 

o 
en 



% 



a 'This author was also a graduate student in the Department of Electrical and Computer Engineering and the 
Department of Physics and Astronomy at Michigan State University, East Lansing, MI 48825 during the preparation 
of this manuscript. 



I. INTRODUCTION 

Langevin dynamics^ is a modeling technique, in which the motion of a set of massive bodies in the 
presence of a bath of smaller solvent particles is directly integrated, while the dynamics of the solvent 
are "averaged out". This approximation may lead to a dramatic reduction in computational cost 
compared to "explicit solvent" methods, since dynamics of the solvent particles no longer need to 
be fully resolved. With Langevin dynamics, the effect of the solvent is reduced to an instantaneous 
drag force and a random delta-correlated force felt by the massive bodies. This framework has 
been applied in a variety of scenarios, including implicit solvents,^ Brownian dynamics, 4 dynamic 
thermostats,— and the stabilization of time integratorsP 

In spite of the numerous successes of conventional Langevin dynamics it has long been rec- 
ognized that there are physically compelling scenarios in which the underlying assumptions break 
down, necessitating a more general treatment P^ To this end, generalized Langevin dynamics (GLD) 
permits the modeling of systems in which the inertial gap separating the massive bodies from the 
smaller solvent particles is reduced. Here the assumptions of an instantaneous drag force and a 
delta-correlated random force become insufficient, leading to the introduction of a temporally non- 
local drag and a random force with non-trivial correlations. GLD has historically been applied to 
numerous problems over the years,™^ with a number of new applications inspiring a resurgence of 
interest, including microrheologyp^"^ biological systemsp*"^ and other situations in which anoma- 
lous diffusion arises.^ 

To facilitate computational exploration of some of these applications, a number of authors have 
developed numerical integration schemes for the generalized Langevin equation (GLE), either in 
isolation or in conjunction with extra terms accounting for external or pairwise forces as would be 
required in a molecular dynamics (MD) simulation. These schemes must deal with a number of 
complications if they are to remain computationally efficient and accurate. 

• The retarded drag force takes the form of a convolution of the velocity with a memory kernel. 
This requires the storage of the velocity history, and the numerical evaluation of a convolution 
at each time step, which can become computationally expensive. 

• The generation of a random force with non-trivial correlations may also require the storage of 
a sequence of random numbers, and some additional computational expense incurred at each 
time step. 

Numerous methods exist that circumvent either one or both of these dimcultiesP^"^ Each has a 



different computational cost, implementation complexity, order of convergence, and specific features, 
e.g., some are restricted to single exponential memory kernels, require linear algebra, etc. To this 
end, it is difficult to distinguish any individual method as being optimal, especially given the broad 
range of applications to which GLD may be applied. 

Motivated by the aforementioned applications and previous work in numerical integrators, we 
have developed a new family of time integration schemes for the GLE in the presence of conservative 
forces, and implemented it in a public domain MD code, LAMMPSP^ Our primary impetus was to 
enable the development of reduced order models for implicit solvents with long-time memory effects. 
Otherwise, to date the GLE has only been solved in the presence of a number of canonically tractable 
external potentials and memory kernelsP 6 | 29 | 3(J | Integration into the LAMMPS framework, provides 
a number of capabilities. LAMMPS includes a broad array of external, bonded, and non-bonded 
potentials, yielding the possibility for the numerical exploration of more complex systems than have 
been previously studied. Finally, LAMMPS provides a highly scalable parallel platform for studying 
N-body dynamics. Consequently, extremely large sample statistics are readily accessible even in 
the case of interacting particles, for which parallelism is an otherwise non-trivial problem. 

In this paper, details germane to both the development of our time integration scheme, as well 
as specifics of its implementation are presented. The time integration scheme to be discussed is 
based upon a two-parameter family of methods specialized to an extended variable formulation of 
the GLE. Some of the salient advantages of our formulation and the final time integration scheme 
are: 

• Generalizability to a wide array of memory kernels representable by a positive Prony series, 
such as power laws. 

• Efficiency afforded by an extended variable formulation that obviates the explicit evaluation 
of time convolutions. 

• Inexpensive treatment of correlated random forces, free of linear algebra, requiring only one 
random force evaluation per extended variable per timestep. 

• Exact conservation of the first and second moments of either the integrated velocity distribu- 
tion or position distribution, for harmonically confined and free particles. 

• Numerically stable approach to the Langevin limit of the GLE. 

• Simplicity of implementation into an existing Verlet MD framework. 



The specialization to Prony-representable memory kernels is worth noting, as there is a growing 
body of literature concerning this form of the GLEP3EQ32] a number of results have been presented 
that establish the mathematical well-posedness of this extended variable GLE including a term 
accounting for smooth conservative forces, as may arise in MD™ These results include proofs of 
ergodicity and exponential convergence to a measure, as well as a discussion of the Langevin limit of 
the GLE in which the parameters of the extended system generate conventional Langevin dynamics. 

One somewhat unique feature of our framework is that we are analyzing the GLE using methods 
from stochastic calculus. In particular, we focus on weak convergence in the construction our 
method, i.e., error in representing the moments of the stationary distribution or the distribution 
itself. The optimal parametrization of our two-parameter family of methods will be defined in terms 
of achieving accuracy with respect to this type of convergence. In particular, the optimal method 
that has been implemented achieves exactness in the first and second moments of the integrated 
velocity distribution for harmonically confined and free particles. Few authors have considered 
this type of analysis for even conventional Langevin integrators, with a notable exception being 
Wang and Skeelj^ who have carried out weak error analysis for a number of integrators used in 
conventional Langevin dynamics. To the best of our knowledge, this is the first time that such a 
weak analysis has been carried out for a GLE integrator. We hope that these considerations will 
contribute to a better understanding of existing and future methods. 

The remainder of this paper is structured as follows: 

• Section III Al introduces the mathematical details of the GLE. 



Section II B presents the extended variable formulation and its benefits. 



Section II C develops the theory associated with integrating the extended variable GLE in 



terms of a two-parameter family of methods. 



Section II D| provides details of the error analysis that establishes the 'optimal' method 



among this family. 



Section III summarizes details of the implementation in LAMMPS. 



Section IV presents a number of results that establish accuracy in numerous limits/scenarios, 



including demonstration of utility in constructing reduced order models. 



II. MATHEMATICAL DETAILS 

A. Statement of the Problem 

The GLE for N p particles moving in ci-dimensions can be written as 

t 
MdV(t) = F c (X(t)) dt- J T(t - s)V(s)dsdt + F r (t)dt, (la) 

o 
dX(t) = V(t)dt, (lb) 

with initial conditions X(0) = X and V(0) = Vq. Here, F c : R. dNp — y R. dNp is a conservative force, 
F r is a random force, M is a diagonal mass matrix, and T is a memory kernel. The solution to this 
stochastic integro-differential equation is a trajectory, which describes the positions X : M + — y M> dNp 
and velocities V : M. + — y M> dNp of the particles as a function of time, t > 0. The second term on 



the right-hand side of Equation [Ta] accounts for the temporally non-local drag force, and the third 
term accounts for the correlated random force. The nature of both forces are characterized by the 
memory kernel, T : 1R + — y R, consistent with the Fluctuation-Dissipation theorem (FDT)P22 The 
FDT states that equilibration to a temperature, T, requires that the two-time correlation of F r (t) 
and r(t) be related as follows: 

(F[(t + s)FJ(t)} = k B TT(s)8 ij , s > 0. (2) 

Here, Sij is the Kronecker delta, and k^ is Boltzmann's constant. 

In the context of an MD simulation, we are interested in solving Equation [I] for both X{t) and 
V(t) at a set of N t uniformly spaced discrete points in time. To this end, we seek to construct a 
solution scheme that is mindful of the following complications: 

• Calculation of the temporally non-local drag force requires a convolution with T(t), and thusly 
the storage of some subset of the time history of V(t). 

• Numerical evaluation of F r (t) requires the generation of a sequence of correlated random 
numbers, as specified in Equation |2j 

• As Equation [T] is a stochastic ODE, we are not concerned with issues of local or global error, 
but rather that the integrated solution converges in distribution. 

To circumvent the first two complications, we work with an extended variable formalis m 1 14 * 35 * 36 ^ in 



which we assume that T(t) is representable as a Prony series: 

N k 

Ck 






t 

Tk 



t > 0. 



(3) 



As will be demonstrated in Section II B , this form of the memory kernel will allow us to map the 
non-Markovian form of the GLE in Equation [T] onto a higher-dimensional Markovian problem with 



dNp, extended variables per particle. The third complication is resolved in Sections II C and II D 
in which a family of integrators is derived, and then 'optimal' parameters are selected based upon 
an error analysis of the moments of the integrated velocity 



B. Extended Variable Formalism 



We introduce the extended variable formalism in two stages. First, we define a set of extended 
variables that allow for an effectively convolution-free reformulation of Equation [T] Then, we 
demonstrate that the non-trivial temporal correlations required of F r (t) can be effected through 
coupling to an auxiliary set of Ornstein-Uhlenbeck (OU) processes. 

We begin by defining the extended variable, Z^it), associated with the fcth Prony mode's action 

on the ith component of X(t) and V(t): 

t 

(t-s 



Zi,k(t) 



Ck 

Tk 



exp 



Tk 



Vi(s)ds 



(4) 



Component-wise, Equation [T] can now be rewritten as: 

N k 

rriidViit) = F t c (X(t)) dt + J2 z i,k(t)dt + F[{t)dt, 

k=l 



(5a) 



dXi(t) = Vi(t)dt. (5b) 

Rather than relying upon the integral form of Equation [4] to update the value of Zi jk {t), we consider 
the total differential of Zi^it) to generate an equation of motion that takes the form of a simple 
SDE: 

dZ itk (t) = --Z hk (t)dt - -Vi(t)dt (6) 

Tk T k 

Now, the system of Equations [5] and |6J can be resolved for X^t), Vi(t), and Z itk (t) without requiring 
the explicit evaluation of a convolution integral. 

Next, we seek a means of constructing random forces that obey the FDT, as in Equation [2] To 
this end, we consider the following SDE: 

dF hk {t) = -—F itk (t)dt + ^y/2k B Tc k dW itk {t) (7) 



Tk 



Tk 



If Wi t k is a standard Wiener process, this SDE defines an Ornstein-Uhlenbeck (OU) process, F itk (t). 
Using established properties of the OU process,^ we can see that F i:k (t) has mean zero and two-time 
correlation: 

(F iik (t + s)F iik (t)) = k B T-exp s, s>0 (8) 

Tfe I T k . 

It is then clear, that the random force in Equation [T] can be rewritten as: 



F[(t) = J2^At) (9) 

fc=i 

Here each individual contribution is generated by a standard OU process, the discrete-time version 
of which is the AR(1) process. While we are still essentially forced to generate a sequence of 
correlated random numbers, mapping onto a set of AR(1) processes has the advantage of requiring 
the retention of but a single prior value in generating each subsequent value. Further, standard 
Gaussian random number generators can be employed. 

Combining both results, the final extended variable GLE can be expressed in terms of the 
composite variable, Si jk (t) = Z^ k (t) + F itk (t): 

N k 

m l dV l {t) = Ft (X(t)) dt + J2 S i,kdt (10a) 

fe=i 

dXiit) = Vi(t)dt (10b) 

dS i>k (t) = --S hk (t)dt - -V t (t)dt + -^2k B Te k dWi, k (t) (10c) 

n Tk Tk 



It is for this system of equations that we will construct a numerical integration scheme in Section II C 
It is worth noting that other authors have rigorously shown that this extended variable form of 
the GLE converges to the Langevin equation in the limit of small r k . 32 Informally, this can be seen 
by multplying the S{ k equation by T k , and taking the limit as r k goes to zero, which results in 



Si, k (t)dt = -c k Vi(t)dt + y/2k B Tc k dW itk (t). 

Inserting this expression into the equation for Vi, we obtain 

(N k \ N k 

^c k \v l (t)dt + ^K^fJ^^dw^it), 
k=l J fe=l 

dXi(t) = Vi(t)dt, 

which is a conventional Langevin equation. We have been careful to preserve this limit in our 
numerical integration scheme, and will explicitly demonstrate this theoretically and numerically. 



C. Numerical Integration of the Extended GLE 



We consider a family of numerical integration schemes for the system in Equation 10 assuming 
a uniform timestep, At. Notation is adopted such that X^nAt) = X™ for hgN. Given the values 
of Xf, V™, and S™ k , we update to the (n + l)th time step using the following splitting method: 



1. Advance Vi by a half step: 

Ay Ay. N h 

+1/2 = m m j, 

2. Advance X{ by a full step: 

X n+1 = X n + AtT /.«+l/2 ^ 

3. Advance <%,& by a full step: 

S^t 1 = OkSb - (1 - B h ) c k V; +1/2 + a k ^2k^fc~ k Bl k (13) 

4. Advance Vi by a half step: 

Ay Ay N h 

vr +l = vr 1/2 + ^-f; (x^) + ^ £ s ?f ( 14 ) 



2m, 2m; 

1 l fc=i 



Here, each B™ k is drawn from an independent Gaussian distribution of mean zero and variance 
unity. The real-valued 6 k and a k can be varied to obtain different methods. For consistency, we 
require that 



6 k = l- — + £>(A£ 2 ), and a k = ^— * + 0(At). 

T k r k 

For the remainder of this article, we restrict our attention three different methods, each of which 
corresponds to a different choice for 9 k and a k . 

• Method 1: Using the Euler-Maruyama scheme to update S^ k is equivalent to using 

At y/M 

U k := 1 and a k := 

Tk T k 

• Method 2: If Vi is held constant, the equation for S^k can be solved exactly. Using this 
approach is equivalent to setting 



'1 - 9 2 ) 
6 k := exp(-At/r fc ) and a k := ^/^— — - 



Method 3: Both methods 1 and 2 are unstable as r k goes to zero. To improve the stability 
when Tfc is small, we consider the following modified version of method 2: 



l-^ 2 



k := exp(-At/r fe ) and a k := \l v ^ k/ 

Note that all three methods satisfy the consistency condition, and are equivalent to the Stormer- 
Verlet-leapfrog method^ when c k = and S^O) = 0. 

D. Error and Stability Analysis 

To help guide our choice of method, we compute the moments of the stationary distribution for a 
one-dimensional harmonic potential (natural frequency u) and a single mode memory kernel (weight 
c and time scale r) . A similar approach has been used for the classical Langevin equation] 33 * 39 ! The 
extended variable GLE for this system converges to a distribution of the form 

p(X, V,S) = -= exp \-(mV 2 + mu 2 X 2 + -S 2 )/2k B T 
Z L c . 

Where Z is the usual normalization constant. From this, we can derive the analytic first and second 
moments 

(V) = 0, (X) = 0, and (S) = 0. 

{XV} = 0, (VS) = 0, and (XS) = 0. 

(V 2 ) = k ^, (X 2 ) = ^4, and (S 2 ) = °^. 

Next, we consider the discrete-time process generated by our numerical integrators, and show 
that the moments of its stationary distribution converge to the analytic ones in At. 

yn+l/2 = yn_ ^ 2 jn + ^gn 

2 2m 

X n+l = X n + At V H+1/2 

S n+i = gs n -(l- 9)cV n+1/2 + ay/2k B TcB n 

yn+l _ yn+l/2 _ ^Ji^n+l , ^ gn+1 

2 2m 

The stationary distribution of this process is defined by the time independence of its first moments 

(X n+1 ) = (X n ), (V n+1 ) = (V n ), and (S n+1 ) = (S n ) 

Enforcing these identities, it can be shown that 

(V n ) = 0, (X n ) = 0, and (S n ) = 0. 



Hence, the first moments are correctly computed by the numerical method for any choice of 9 and 
a. Computing the second moments, we obtain 

(*.y> = o, (vs-Ho, (*-*-> = 2cAt(1 _ e) ^%% {i _ (uAm ■ 

m 2 _ k B T Ato 2 2 _ k B T Ato 2 (2cAt(l - #) 2 - 4m(l - # 2 )) 



<(V»)') = -?_ 7 -J^-, ((^T) 



m (1-fl) 2 ' xv y/ mw 2 (l-#) 2 (2cAt(l-#) 2 -m(l-# 2 )(4-(u;At) 2 ))' 
and ((S»)V ek » T Wm < 4 -^ 



r m(l-9 2 )(A-(cuAt) 2 )-2cAt(l-9) 2 
From this analysis, we conclude that we obtain the correct second moment for V for any method 
with Ato 2 = (1 — 9) 2 . Now, applying the particular values of 9 and a, and expanding in powers of 
At, we obtain the following. 



Method 1: 



<(n 2 > = — , ((x n f) = ^ (i + ^^ ) + o(At^ 



m '' ' ' mco 2 \ 4 

(X n S n ) = At ] CkBT +0(At 3 ), and ((S n ) 2 ) = — (l + ^ )+ 0(At 2 ) 
Amr t \ 2t J 



Method 2: 



rn ^ k B T / At 2 \ , „,. 4 . „_. 2 . k B T / ( (l + 3(a;r) 2 )At 2 . , 



At 2 ck B T , ml . 4 2 ck B T ^ , cAt 2 \ 



(X^) = ^^ + (9(At 4 ), and <(^) 2 > = -^ ^1 + —j + O(Af) 
• Method 3: 

(( v)»> = *£, <m'> = ig (i + <^) + (Af), 

<X» 5 ») = ^M + (A^), and ((fl .)» ) = £^:( 1+ (gg^g ) +0( Af) 

For methods 1 and 3, we obtain the exact variance for V, independent of At, since they both satisfy 
Ata 2 = (1 — 9) 2 . For method 2, the error in the variance of V is second-order in At. All three 
methods overestimate the variance of X, with an error which is second-order in At. The error in 
the variance of S is first-order for method 1, and second-order for methods 2 and 3. 

It is possible to choose 9 and a to obtain the exact variance for X, but this would require using 
a different value for 9 and a for each value of u. This is not useful in our framework, since the 
method is applied to problems with general nonlinear interaction forces. 

10 



We would like our numerical method to be stable for a wide range of values for rfc. As we 



mentioned in Section II B the GLE converges to the conventional Langevin equation as r k goes 
to zero, and we would like our numerical method to have a similar property. For fixed At, both 
methods 1 and 2 are unstable (a k is unbounded) as r k goes to zero. However, method 3 does not 
suffer from the same problem, with 9 k converging to zero and a k bounded. 

From this analysis, we conclude that method 3 is the best choice for implementation. Prior to 
providing implementation details, however, we briefly consider a simple multistage extension that 
can capture the exact first and second moments of position and velocity simultaneously at the 
expense of introducing a numerical correlation between them. 

E. Multistage Splitting 

Inspired by the work of Leimkuhler and Matthewsp^ we consider a generalization of the splitting 



method considered in Section II C in which the position and velocity updates are further split, 

2m,- l y ~° ' 2m,- 



Vr l/4 = V? + ^-Ft (X") + (1 - $£- £ s- k 



k=l 



n+1/2 = X n + A^n+1/4 



vr /2 =vr i,4 +^f: s i k 



2m,- 



fc=l 



S?t l = OkSb, - (1 - 6 k ) c k V? +1/2 + a k y/2k B Tc k Bl 



n+3/4 n+1/2 ^|^n +1 

1 k=i 
X n+l =- x n+1 ^ 2 + ^T/ n + 3 / 4 

Ay- Ay. N * 

+1 = +3 /4 M n+1 _ m y* +1 

1 l 2rrii v ' 2m,- *-^ hk 

1 l fe=i 

In the special case that £ = 0, we have V™ = V™ =Vi , the two updates of X can be 
combined, and we recover our original splitting method. 



Repeating the analysis in Section II D for the harmonic oscillator with a single memory term, 
we find 

k B T Ato 2 (cAt(l - 6)\A - (coAtfe) - 2m(l - fl 2 )(4 - (cAt) 2 Q) 
U J ; " mw 2 (1 - #) 2 (cAt(l - 0) 2 (4 - (wA£) 2 - 2m(l - # 2 )(4 - (wAt) 2 )) ' 

11 



mcoHi-er {XV) - 


= 


ckeT 4ra 2 m 





/^ncynv = 4At 2 a 2 (l-Qck B T 

{ } cA£(l-#) 2 (4-(wAt) 2 0-2m(l-fl 2 )(4-(cc;A£) 2 )' 

anH /r^ 2 \ = ^I 4ra 2 m(4-(g;At) 2 ) 

U J ; " r 2m(l - # 2 )(4 - (wAt) 2 ) - cAt(l - #) 2 (4 - (wAt) 2 £) ' 

In the special case that £ = 1, we can simplify these expressions to obtain 
,(V^\ k B T Ato 2 (4-(Ato;) 2 ) 

<*»*»> = 0, <V»S»> = 0, and <(^) 2 > r 2w(1 _ , 2) _ cAi(1 _ , )2 - 

From this analysis, we conclude that we obtain the correct second moment for X for any method 
with £ = 1 and Ata 2 = (1 — 6) 2 . To guarantee the correct second moment for V, the choice of 6 
and a becomes u dependent. As was discussed in the previous section, the parameters prescribed 
by methods 1 and 3 using the original splitting have a similar behavior but with the roles of X and 
V reversed. 

It is then tempting to formulate a method that exactly preserves the second moments of both 
X and V at the same time. It turns out that this is possible by simply shifting where we observe 
V, using either V™ or V™ in the multistage splitting method above. For example, consider 
the following asymmetric method, with £ = 1, 

xr 1 ' 2 = *r + ^vr 

At Nk 



vr i/2 =vr+^-j:^ 



2rrij 
1 fe=i 



«u +1 = «*s« - (i - h) c t v; M ' 2 + arv /2k B T Ci a;; t 



r n+3/4 _ T/ n+l/2 _ At 

2m,- 

1 k=i 

X n+1 — Y n+1 ^ 2 + ^ T/ n + 3 / 4 

n+1 = n+3/4 + At . n+1 . 

rrii 



For this method we obtain, 

/n/^2\ _ k -B T Ata 2 2 _ k B T Ata 2 _ -At 2 a 2 k B T 

((y } } " ^T(T^P' ((X } } " m^(T^p' {X V ' ~~- 2m(l - tf) 2 ' 

(X n S™) = 0, (F n S") = 0, and ({S n ) 2 ) 



,,,.,, ckgT 4ra 2 m 



r 2m(l - # 2 ) - cAt(l 
12 



If we use a method with Ata 2 = (1 — 6) 2 , we find that we obtain the exact moments for X and 
V, but we have introduced an O(At) correlation between X and V. This is in contrast to the 
symmetric methods where this correlation is identically zero. 

III. IMPLEMENTATION DETAILS 



Method 3, as detailed in Section II D has been implemented in the LAMMPS software package. It 
can be applied in conjunction with all conservative force fields supported by LAMMPS. There are a 
number of details of our implementation worth remarking on concerning random number generation, 
initial conditions on the extended variables, and the conservation of total linear momentum. 

The numerical integration scheme requires the generation of Gaussian random numbers, by way 



of B^ k in Equation 13 By default, all random numbers are drawn from a uniform distribution with 
the same mean and variance as the formally required Gaussian distribution. This distribution is 
chosen to avoid the generation of numbers that are arbitrarily large, or more accurately, arbitrarily 
close to the floating point limit. The generation of such large numbers may lead to rare motions 
that result in the loss of atoms from a periodic simulation box, even at low temperatures. Atom 
loss occurs if, within a single time step, the change in one or more of an atom's position coordinates 
is updated to a value that results in it being placed outside of the simulation box after periodic 
boundary conditions are applied. A uniform distribution can be used to guarantee that this will 
not happen for a given temperature and time step. However, for the sake of mathematical rigor, 
the option remains at compile-time to enable the use of the proper Gaussian distribution with the 
caveat that such spurious motions may occur. Should the use of this random number generator 
produce a trajectory in which atom loss occurs, a simple practical correction may be to use a 
different seed and/or a different time step in a subsequent simulation. It is worth noting that the 
choice of a uniform random number distribution has been rigorously justified by Diinweg and PauP^ 
for a number of canonical random processes, including one described by a conventional Langevin 
equation. We anticipate that a similar result may hold for the extended variable GLE presented in 
this manuscript. 

With respect to the initialization of the extended variables, it is frequently the case in MD that 
initial conditions are drawn from the equilibrium distribution at some initial temperature. Details 
of the equilibrium distribution for the extended system are presented in Section 2 of an article by 
Ottobre and PavliotisP^ In our implementation, we provide the option to initialize the extended 
variables either based upon this distribution, or with zero initial conditions (i.e., the extended 

13 



system at zero temperature). As it is typically more relevant for MD simulations, the former is 
enabled by default and used in the generation of the results in this paper. 

Conservation of the total linear momentum of a system is frequently a desirable feature for 
MD trajectories. For deterministic forces, this can be guaranteed to high precision through the 
subtraction of the velocity of the center of mass from all particles at a single time step. In the 
presence of random forces such as those arising in GLD, a similar adjustment must be made at each 
time step to prevent the center of mass from undergoing a random acceleration. While it is not 
enabled by default, our implementation provides such a mechanism that can be activated. When 
active, the average of the forces acting on all extended variables is subtracted from each individual 
extended variable at each time step. While this is a computationally inexpensive adjustment, it 
may not be essential for all simulations. 

IV. RESULTS 

Throughout this section, results will be presented, primarily in terms of the integrated velocity 
autocorrelation function (VAF). This quantity is calculated using "block averaging" and "subsam- 
pling" of the integrated trajectories for computational convenience.^ Error bars are derived from 
the standard deviation associated with a set of independently generated trajectories. 

We begin by presenting results that validate our time integration scheme for a GLE in the absence 
of a conservative force with a single mode Prony series kernel. In this case, the GLE is analytically 
soluble.^ We consider the normalized VAF as a metric for comparison. For a kernel of the form 



r(£) = - exp 



T 



£>0, (15) 



the normalized VAF takes the following form: 

(V(t)V(0)} Jexp[^](cos(fit) + ^sin(fi£)) for ^ 

"«1«p[-a(iH) *»«=»■ a= ^^ r - (16) 

Making an analogy with the canonical damped harmonic oscillator, we consider three scenarios, i. 
underdamped (real Q), ii. critically damped (Q = 0), and iii. overdamped (imaginary fl). In Figure 
[TJ we demonstrate that we can recover all three regimes using our integrator. 

To validate our integration scheme in the presence of a conservative force, we next consider GLD 
with an external potential. The analytic solution of the GLE with a power law memory kernel 
and a harmonic confining potential has been derived by other authors.^ In the cited work, the 
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FIG. 1. Normalized VAF for a single Prony series mode in the i.) underdamped (c = 1, r = 1), ii.) 
critically damped (c = 0.5, r = 0.5), and iii.) overdamped limits (c = 0.25, r = 0.25). A time step of 
At = 0.01 is used for all runs, and error bars are drawn based upon a sample of 10,000 walkers over 10 
independent runs. 

GLE is solved in the Laplace domain, yielding correlation functions given in terms of a series of 
Mittag-Leffler functions and their derivatives. Here, we apply our integration scheme to a Prony 
series fit of the power law kernel and demonstrate that our results are in good agreement with the 
exact result over a finite time interval. The GLE that we intend to model has the form: 

t 
dV{t) = -u%X{t)dt - / 7A (t - s)- x V{s)dsdt + M~ x F r {t)dt, 0<A<1 (17) 

j r(i - A) 

o 

We begin by constructing a Prony series representation of the memory kernel. As it exhibits a 
power law decay in the Laplace domain, as well, a Prony series representation will have strong 
contributions from modes decaying over a continuum of time scales. To this end, rather than 
relying on a non-linear fitting procedure to choose values for r&, we assume logarithmically spaced 
values from At/10 to 10N t At. By assuming the form of each exponential, the Prony series fit 
reduces to a simple linear least squares problem, that we solve using uniformly spaced data over an 
interval that is two decades longer than the actual simulation. In Figure [2j the Prony series fit of 
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the memory kernel for A = 0.5 is compared to the its exact value. 



Exact 
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FIG. 2. A Prony series fit of the power law memory kernel in Equation 17 for 7^ = 1, A = 0.5 with an 
increasing number of modes. 

Here, the maximum relative error is ~ 10%, while it is ~ 1% for Nk = 8. In Figure |3j the 
normalized VAF computed via numerical integration of the extended GLE with a variable number 
of modes is shown compared to the exact result for some of the parameters utilized in an article by 
Desposito and VinalesP^ It seems evident that the accuracy of the integrated velocity distribution 
improves relative to the exact velocity distribution as the number of terms in the Prony series fit 
increases. This is quantified in Figure |4j in which the pointwise absolute error in the integrated 
VAF is illustrated. 

Next, we demonstrate that our implementation is robust in certain limits of the GLE - particu- 
larly the Langevin and zero coefficient limits. As our implementation is available in a public domain 
code with a large user base, developing a numerical method that is robust to a wide array of inputs 
is essential. For the zero coefficient limit we consider a harmonically confined particle experiencing 
a single mode Prony series memory kernel of the following form: 



dV{t) =- -uj 2 Q X{t)dt - J - ( 



-(t-s)M 



V(s)dsdt + M- L F r (t)dt, < A < 1 



(18) 



The initial conditions on X and V are drawn from a thermal distribution at T = 1. In the limit 
that c — > 0, we expect the integrated normalized VAF to approach that of a set of deterministic 
harmonic oscillators. The initial conditions on X and V are drawn from a thermal distribution at 
T — 1, giving the integrated result some variance, even in the Newtonian/deterministic limit. In 
Figure |5j that this limit is smoothly and stably approached is illustrated. 

In the exact c = limit, oscillations in the VAF occur at the natural frequency of the confining 
potential, Uq — 1, whereas for non-zero c, these oscillations are damped in proportion to c, as one 
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FIG. 3. Comparison of numerically integrated results to the exact solution of Equation 17 for 7^ = 1, 
A = 0.5, and uo = 1.4. A time step of At = 0.01 is used, and error bars are drawn based upon a sample of 
10,000 walkers over 10 independent runs with At = 0.01. 




FIG. 4. Pointwise absolute error in the normalized velocity autocorrelation function for the same conditions 
as Figure [3j Error is computed with respect to the mean of VAFs computed from 10 independent runs. 



may expect on the basis of intuition. The Langevin limit is illustrated next. To this end, we utilize 
the same single mode Prony series memory kernel, but remove the confining potential (i.e., u>o = 0). 
We have done so to ensure that the Langevin limit yields an Ornstein-Uhlenbeck process. In Figure 
[6j we illustrate that this limit is also smoothly and stably approached. 
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FIG. 5. Illustration of the smooth/stable manner in which the proposed integrator approaches the c = 
(Newtonian) limit from c = 1 for a single Prony mode with r = 1. The particle is harmonically confined 
with ooq = 1, and the expected period of oscillation is restored for c = 0. 




FIG. 6. Illustration of the smooth/stable manner in which the proposed integrator approaches the r = 
(Langevin) limit from r = 1 for a single Prony mode. The particle is subject to no conservative forces, 
so the resultant dynamics correspond to an Ornstein-Uhlenbeck process. The 'exact' Langevin limit was 
integrated using fix langevin in LAMMPS. 



The result for the Langevin limit itself was integrated using the existing fix langevin command 
in LAMMPS. For r = 10°, the GLE yields results that differ from the Langevin limit as one might 
expect on the basis of the previous results. However, even for r = lCT 1 , the GLE yields results that 
are close to that of the Langevin limit, as the resultant numerical method retains some GLE-like 
behavior. However, for r = 1CF 8 and At = 0.01, the parameter 9 = exp [— 10 6 ], is below machine 



epsilon. As mentioned in Section |HD[ this results in a method that is indistinguishable from 
conventional Langevin dynamics, and consequently, the integrated dynamics are indistinguishable. 
This result demonstrates that the method is robust, even if the user 'accidentally' crosses into a 
Langevin-like regime. We again emphasize that methods 1 and 2, considered earlier, will not stably 
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approach this limit. 

To demonstrate capability in reduced order modeling, we next invert the VAF associated with a 
set of trajectories generated in the presence of a conservative force, to construct a memory kernel 
that will reproduce the same dynamics without such a force. The effect of the force is to be included 
in the GLE memory kernel. The starting point for this procedure is the observation that the VAF 
and memory kernel have the following relationship in the Laplace domain^ 

k B T 



(V(a)V(O)) 



(19) 



ma + r(cr) 

Given the Laplace domain VAF, we can solve for the Laplace domain memory kernel algebraically. 
Applying an inverse Laplace transform then yields a memory kernel that will deliver the same VAF: 

k B T 



r(t) = £- 1 {- 



ma}. 



(20) 



(V(a)V(O)) 

This memory kernel can then be fit to a Prony series, and utilized in the numerical integration 
scheme being developed to generate trajectories that give rise to the same velocity distribution, but 
without the explicit inclusion of a conservative force. In Figure [7| we demonstrate that we can take 



the exact VAF generated by Equation [TTJ and reproduce its dynamics with a force-free GLE with 
this inversion procedure. 




FIG. 7. Illustration of the result of the kernel fitting procedure described in Section IV Note that while the 
exact VAF evolves under the influence of a harmonic confining potential, the integrated VAF is spawned 
by a GLE free of any conservative forces. The effect of the force is included in the GLE memory kernel 
itself. 



The parameters utilized are ujq = 0.5 and A = 0.9. The VAF has been fit to a three term Prony 



series, which is then Laplace transformed and inserted into Equation 20 The resultant time domain 
memory kernel is then also fit to a 3 term memory kernel, one term of which has been determined 
to be a conventional Langevindike term. The reconstructed memory kernel is then applied to a 
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conservative force-free GLE, and the dynamics are integrated as with all other results. This result 
is especially interesting as it demonstrates capability for effecting a confining potential, on a finite 
timescale, exclusively through the GLE drag/random forces. This same fitting procedure could be 
used to remove inter-particle interactions, albeit with some loss of dynamical information, and will 
be discussed in more detail in future work. 

It is important to note that while the dynamics of a confined particle are reproduced in the 
absence of an explicit confining force, this is only the case over the interval for which the mem- 
ory kernel is reconstructed. Outside of this interval, the asymptotic behavior of the Prony series 
memory kernel will give rise to unbounded diffusive motion. While the expected discrepancies in 
the dynamics are difficult to notice in the VAF, as both the analytic and reconstructed quantities 
decay to zero, it is very evident in the mean-squared displacement (MSD). Here the asymptotic 
exponential behavior of the Prony series memory kernel necessarily leads to asymptotic diffusive 
motion signaled by a linear MSD. In contrast, the analytic memory kernel with the confining force 
will generate an MSD that remains bounded. This is illustrated in Figure [8} Here, the fit VAF is 
integrated to yield the MSD for the Prony series model, and compared with the exact MSD over 
both the interval of the fit, and beyond. This example illustrates the importance of choosing an 
appropriate interval for fitting memory kernels to achieve the appropriate limiting behavior in one's 
dynamics. 
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FIG. 8. Comparison of the exact normalized MSD and the MSD generated from integrating the VAF fit. 
The exact MSD is normalized such that its asymptotic limit is 1, and the same normalization is applied to 
the MSD integrated from the fit. The grey line indicates the upper bound of the interval over which the 
VAF fit was constructed. 
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V. CONCLUSIONS 

A family of numerical integration schemes for GLD have been presented. These schemes are 
based upon an extended variable formulation of the GLE in which the memory kernel is rendered 
as a positive Prony series. In certain limits, it can be shown that a specific instance of this family of 
integrators exactly conserves the first and second moments of the integrated velocity distribution, 
and stably approaches the Langevin limit. To this end, we identify this parametrization as optimal, 
and have implemented it in the MD code, LAMMPS. Numerical experiments indicate that this 
implementation is robust for a number of canonical problems, as well as certain pathologies in the 
memory kernel. An exemplary application to reduced order modeling illustrates potential uses of 
this module for MD practitioners. Future work will further develop the VAF fitting procedure in 
the context of statistical inference methods, and present extensions of the numerical integrator to 
mixed sign and complex memory kernels. 
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